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When uncontrollable resources fluctuate, Optimum Power Flow (OPF), routinely used by the 
electric power industry to re-dispatch hourly controllable generation (coal, gas and hydro plants) 
over control areas of transmission networks, can result in grid instability, and, potentially, cascading 
outages. This risk arises because OPF dispatch is computed without awareness of major uncer- 
tainty, in particular fluctuations in renewable output. As a result, grid operation under OPF with 
renewable variability can lead to frequent conditions where power line flow ratings are significantly 
exceeded. Such a condition, which is borne by simulations of real grids, would likely resulting in 
automatic line tripping to protect lines from thermal stress, a risky and undesirable outcome which 
compromises stability. Smart grid goals include a commitment to large penetration of highly fluc- 
tuating renewables, thus calling to reconsider current practices, in particular the use of standard 
OPF. Our Chance Constrained (CC) OPF corrects the problem and mitigates dangerous renewable 
fluctuations with minimal changes in the current operational procedure. Assuming availability of 
a reliable wind forecast parameterizing the distribution function of the uncertain generation, our 
CC-OPF satisfies all the constraints with high probability while simultaneously minimizing the cost 
of economic re-dispatch. CC-OPF allows efficient implementation, e.g. solving a typical instance 
over the 2746-bus Polish network in 20s on a standard laptop. 

The power grid can be considered one of the greatest engineering achievements of the 20th century, responsible for 
the economic well-being, social development, and resulting political stability of billions of people around the globe. 
The grid is able to deliver on these goals with only occasional disruptions through significant control sophistication 
and careful long-term planning. 

Nevertheless, the grid is under growing stress and the premise of secure electrical power delivered anywhere and at 
any time may become less certain. Even though utilities have massively invested in infrastructure, grid failures, in 
the form of large-scale power outages, occur unpredictably and with increasing frequency. In general, automatic grid 
control and regulatory statutes achieve robustness of operation as conditions display normal fluctuations, in particular 
approximately predicted inter-day trends in demand, or even unexpected single points of failure, such as the failure 
of a generator or tripping of a single line. However, larger, unexpected disturbances can prove quite difficult to 
overcome. This difficulty can be explained by the fact that automatic controls found in the grid are largely of an 
engineering nature (i.e. the flywheel-directed generator response used to handle short-term demand changes locally) 
and are largely not of a data-driven, algorithmic and distributed nature. Instead, should an unusual condition arise, 
current grid operation relies on human input. Additionally, only some real-time data is actually used by the grid to 
respond to evolving conditions. 

All engineering fields can be expected to change as computing becomes ever more enmeshed into operations and 
massive amounts of real-time data become available. In the case of the grid this change amounts to a challenge; 
namely how to migrate to a more algorithmic-driven grid in a cost-effective manner that is also seamless and gradual 
so as not to prove excessively disruptive - because it would be impossible to rebuild the grid from scratch. One of 
the benefits of the migration, in particular, concerns the effective integration of renewables into grids. This issue is 
critical because large-scale introduction of renewables as a generation source brings with it the risk of large, random 
variability - a condition that the current grid was not developed to accommodate. 

This issue becomes clear when we consider how the grid sets generator output in "real time" . This is typically 
performed at the start of every fifteen-minute (to an hour) period, or time window, using fixed estimates for conditions 
during the period. More precisely, generators are dispatched so as to balance load (demand) and generator output 
at minimum cost, while adhering to operating limitations of the generators and transmission lines; estimates of the 
typical loads for the upcoming time window are employed in this computation. This computational scheme, called 
Optimum Power Flow (OPF) or economic dispatch, can fail, dramatically, when renewables are part of the generation 
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FIG. 1: Bonneville Power Administration [l[ shown in outline under 9% wind penetration, where green dots mark actual 
wind farms. We set standard deviation to be 0.3 of the mean for each wind source. Our CC-OPF (with 1% of overload set 
as allowable) resolved the case successfully (no overloads), while the standard OPF showed 8 overloaded lines, all marked in 
color. Lines shown orange are at 4% chance of overload. There are two dark red lines which are at 50% of the overload while 
other (dark orange) lines show values of overload around 10%. 

mix and (exogenous) fluctuations in renewable output become large. By "failure" we mean, in particular, instances 
where a combination of generator and renewable outputs conspire to produce power flows that significantly exceed 
power line ratings. When a line's rating is exceeded, the likelihood grows that the line will become tripped (be 
taken uncontrollably out of service) thus compromising integrity and stability of the grid. If several key lines become 
tripped a grid would very likely become unstable and possibly experience a cascading failure, with large losses in 
serviced demand. This is not an idle assumption, since firm commitments to major renewable penetration are in 
place throughout the world. For example, 20% renewable penetration by 2030 is a decree in the U.S. Q, and similar 
plans are to be implemented in Europe, see e.g. discussions in ■ At the same time, operational margins (between 
typical power flows and line ratings) are decreasing and expected to decrease. 

A possible failure scenario is illustrated in Fig. Q] using as example the U.S. Pacific Northwest regional grid data 
(2866 lines, 2209 buses, 176 generators and 18 wind sources), where lines highlighted in red are jeopardized (flow 
becomes too high) with unacceptably high probability by fluctuating wind resources positioned along the Columbia 
river basin (green dots marking existing wind farms). We propose a solution that requires, as the only additional 
investment, accurate wind forecasts; but no change in machinery or significant operational procedures. Instead, we 
propose an intelligent way to modify the optimization approach so as to mitigate risk; the approach is implementablc 
as an efficient algorithm that solves large-scale realistic examples in a matter of seconds, and thus is only slightly 
slower than standard economic dispatch methods. 

Maintaining line flows within their prescribed limits arises as a paramount operational criterion toward grid stability. 
In the context of incorporating renewables into generation, a challenge emerges because a nominally safe way of 
operating a grid may become unsafe - should an unpredictable (but persistent) change in renewable output occur, the 
resulting power flows may cause a line to persistently exceed its rating. It is natural to assess the risk of such an event 
in terms of probabilities, because of the non-deterministic behavior of e.g. wind; thus in our proposed operational 
solution we will rely on techniques involving both mathematical optimization and risk analysis. 

When considering a system under stochastic risk, an extremely large variety of events that could pose danger 
might emerge. Recent works 0-0] suggest that focusing on instantons, or most-likely (dangerous) events, provides a 
practicable route to risk control and assessment. However, there may be far too many comparably probable instantons, 
and furthermore, identifying such events does not answer the question of what to do about them. In other words, we 
need a computationally efficient methodology that not only identifies dangerous, relatively probable events, but also 
mitigates them. 

This paper suggests a new approach for handling the two challenges, that is to say, searching for the most proba- 
ble realizations of line overloads under renewable generation, and correcting such situations through control actions, 
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simultaneously and efficiently in one step. Our approach relics on methodologies recently developed in the optimiza- 
tion literature and known under the name of "Chance-Constrained" (CC) optimization (9|. Broadly speaking, CC 
optimization problems are optimization problems involving stochastic quantities, where constraints state that the 
probability of a certain random event is kept smaller than a target value. 

To address these goals, we propose an enhancement of the standard OPF to be used in the economic dispatch 
of the controllable generators. We model each bus that houses a power source subject to randomness to include a 
random power injection, and reformulate the standard OPF in order to account for this uncertainty. The formulation 
minimizes the average cost of generation over the random power injections, while specifying a mechanism by which 
(standard, i.e. controllable) generators compensate in real-time for renewable power fluctuations; at the same time 
guaranteeing low probability that any line will exceed its rating. This last constraint is naturally formulated as a 
chance constraint - we term out approach Chance-Constrained OPF, or CC-OPF. 

This paper is organized as follows. In "Models" we first describe the various mathematical models used to describe 
how the grid operates, as well as our proposed methodology. We then present, in "Experiments" a number of examples 
to demonstrate the speed and usefulness of our approach. Finally, in "Methods" wc provide in-depth technical details 
on our approach. Further details, and proofs, arc provided in the "Supporting Information" (SI). 

MODELS 

Transmission Grids: Controls and Limits 

The power systems we consider in this paper are transmission grids which operate at high voltages so as to convey 
power economically, with minimal losses, over large distances. This is to be contrasted with distribution systems; 
typically residential, lower voltage grids used to provide power to individual consumers. From the point of view of 
wind-power generation, smooth operation of transmission systems is key since reliable wind sources are frequently 
located far away from consumption. 

Transmission systems balance consumption/load and generation using a complex strategy that spans three different 
time scales (see e.g. (loj). At any point in time, generators produce power at a previously computed base level. Power 
is generated (and transmitted) in the Alternating Current (AC) form. An essential ingredient toward stability of the 
overall grid is that all generators operate at a common frequency. In real time, changes in loads are registered at 
generators through (opposite) changes in frequency. A good example is that where there is an overall load increase. In 
that case generators will marginally slow down - frequency will start to drop. Then the so-called primary frequency 
control, normally implemented on gas and hydro plans with so-called "governor" capability will react so as to stop 
frequency drift (large coal and nuclear units are normally kept on a time constant output). This is achieved by 
having each responding generator convey more power to the system, proportionally to the frequency change sensed. 
(In North America the proportionality coefficient is normally set to 5% of the generator capacity for 0.5Hz deviation 
from the nominal frequency of 60Hz.) This reaction is swift and local, leading to stabilization of frequency across the 
system, however not necessarily at the nominal 60Hz value. The task of the secondary, or Automatic Gain Control 
(AGC), is to adjust the common frequency mismatch and thus to restore the overall balance between generation and 
consumption, typically in a matter of minutes. Only some of the generators in a local area may be involved in this 
step. The final component in the strategy is the tertiary level of control, executed via the OPF algorithm, typically 
run as frequently as every fifteen minutes (to one hour), and using estimates for loads during the next time window, 
where base (controllable) generator outputs are reset. This is not an automatic step in the sense that a computation 
is performed to set these generator levels; the computation takes into account not only load levels but also other 
parameters of importance, such as line transmission levels. Tertiary control computation, which is in the center of 
this paper, thus represents the shortest time scale where actual off-line and network wide (in contrast to automatic 
primary and secondary controls of frequency) optimal computations are employed. The three levels are not the only 
control actions used to operate a transmission system. Advancing further in the time scale, OPF is followed by the 
so-called Unit Commitment (UC) computation, which schedules the switching on and off of large generation units on 
the scale of hours or even days. 

A critical design consideration at each of the three control levels is that of maintaining "stability" of the grid. The 
most important ingredient toward stable operation is synchrony - ultimately, all the generators of the network should 
stabilize thus locking, after a perturbation followed by a seconds-short transient, at the same frequency. Failure to 
do so not only proves inefficient but, worse, it threatens the integrity of the grid, ultimately forcing generators to 
shut down for protective reasons - thus, potentially, causing a large, sudden change in power flow patterns (which 
may exceed equipment limits, see below) and possibly also an unrecoverable generation shortage. A second stability 
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goal is that of maintaining large voltages. This is conducive to efficiency; lower voltage levels cause as a byproduct 
more generation (to meet the loads) and larger current values. Not only is this combination inefficient, but in an 
extreme case it may make impossible to meet existing loads (so-called "voltage collapse" is a manifestation of this 
problem). The third stability goal, from an operational perspective, is that of maintaining (line) power flows within 
established bounds. In long transmission lines, a large flow value will cause excessive voltage drop (an undesirable 
outcome as discussed). On a comparatively shorter line, an excessively large power flow across the line will increase 
the line temperature to the point that the line sags, and potentially arcs or trips due to a physical contact. For each 
line there is a given parameter, the line rating (or limit) which upper bounds flow level during satisfactory operation. 

Of the three "stability" criteria described above, the first two (maintaining synchrony and voltage) are a concern 
only in a truly nonlinear regime which under normal circumstances occur rarely. Thus, we focus on the third - 
observing line limits. 



OPF — Standard Generation Dispatch (tertiary control) 



OPF is a key underlying algorithm of power engineering, see the review in [llj and e.g. [10|, [L2J. The task of OPF, 
usually executed off-line at periodic intervals, is to re-dispatch generation over a control area of the transmission grid, 
for example over the Bonneville Power Administration (BPA) grid shown in Fig. [TJ In outline the standard OPF can 
be stated as the following constrained optimization problem: 



OPF 

PF 
TL 
GB 



min c(p), s.t. (1) 
p 

Power Flow Equations 
Line flow (capacity) Limits 
Generation Bounds 



Here, p = (pi\i S S) is the vector of controllable/re-dispatchable generation available at the subset S of nodes of the 
full set of grid nodes, V; S C V. The above problem is endowed by three types of constraints: Power flow (PF), Line 
Limit (TL) and Generation Bound (GB) constraints. (PF) consists of AC Power Flow equations (AC-PF) which are 
simply Kirchoff's circuit laws stated in terms of power flows and potentials (voltages). Here, for each node i € V its 
voltage Ui is defined as w i e j6 ' i , where Vi and 8i are the voltage magnitude and phase angle at node i. See e.g. 12 |- 



The power flow equations are quadratic and thus can constitute an obstacle to solvability of OPF (from a technical 
standpoint, they give rise to nonconvexities) . In transmission system analysis a linearized version of the AC equations 
is commonly used, the so-called "DC-approximation" . In this approximation (a) all voltages are assumed fixed and 
re-scaled to unity; (b) phase differences between neighboring nodes are assumed small, V{i,j} 6 £ : \9i — 9j\ <C I, 
where £ stands for the set of the grid edges, or lines; (c) thermal losses arc ignored (reactance dominates resistance for 
all lines). Then, the power flow over line {i,j}, with line susceptance ftij (= (3ji) is related linearly to the respective 
phase difference, fij = Pij(9i — 0j), while the balance of power can be stated in the following matrix form 

PF : B6 =p-d, where B = (Bij\i,j e V), (2) 

[ 0, otherwise 

where d is the vector of the exogenous (uncontrollable) demands at each node (possibly equal to zero at some nodes). 
We also model an uncontrollable demand as a negative load. The linearized approximation is often considered around 
a stationary solution of the full AC-PF system (stationary operational point); thus elements of the susceptance matrix 
B account for renormalization due to the base case solution. The GB and TL constraints in problem (JTJ) are 

GB : VieS: pT' m < Pt < P? ax , (4) 
TL: V{i,j}G£: |/y| < (5) 

where p mm i p max are lower and upper generation bounds. f max represents the line limit (typically a thermal limit), 
which is assumed to be strictly enforced in constraint (TL). This conservative condition will be relaxed in the following. 
Finally, the objective c(p) to minimize in problem ([TJ) is a sum of convex quadratic functions of the components of 
p (fixed price curve per generator). In summary, problem ([TJ) is a convex optimization problem solved for a fixed 
vector of demands, d. In practice, however, demand will fluctuate around d; generators then respond by adjusting 
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their output proportionally to the overall fluctuation as discussed above in relation to frequency control. When some 
of the generation is due to renewables, standard OPF would model their output as constant (in the time window of 
interest), and would manage their fluctuation by having controllable generators adjust their output in the same way 
used to handle demand variations. 

Using modern optimization tools (p} is an easily solved problem. This scheme works well in current practice, as 
demands do not substantially fluctuate on the time scale for which OPF applies. Thus the standard practice of 
solving ([TJ in the feasibility domain defined by Eq. ([2j) , (J3j) , (J4j) , ([5]) , using demand forecasts based on historical data 
(and ignoring fluctuations) has produced a very reliable result - generation re-dispatch covering a span of fifteen 
minutes to an hour, depending on the system. 



Chance constrained OPF: motivation 



The separation of generation control into the hierarchy of primary, secondary and OPF has worked in the past 
because of the slow time scales of change in uncontrolled resources (mainly loads). That is to say, frequency control 
and load changes were well-separated. Clearly, an error in the forecast or an under-estimation of possible d for the next 
-e.g., fifteen minute- period may lead to an operational problem in standard OPF; see e.g. the discussions in [l3, 14 1. 
This was not considered a significant handicap until recently, however, simply because line trips due to overloading as 
a result of OPF-directed generator dispatch were (and still are) rare. The projected increase of renewable penetration 
in the future, accompanied by the decreasing gap between normal operation and limits set by line capacities, will 
make these overload events more frequent and generally increase risk (see 0] ) . One would also suspect that the rare 
overload event in the grid of today is due to setting up the TL limits too conservatively. In general, lowering the 
TL limits succeeds in preventing overloads, but it also forces excessively conservative choices of t he g eneration re- 
dispatch, potentially causing extreme volatility of the electricity markets. (See e.g. the discussion in [15| on abnormal 
price fluctuations in ERCOT and New Zealand markets, which are both heavily reliant on renewables.) CC-OPF, 
introduced next, is less conservative (it is probabilistic) and also offers an exact and efficient algorithm for balancing 
the cost of operation with the risk of overload. We will introduce uncertain power sources into the OPF. This will 
change the optimization problem from deterministic to probabilistic; we will seek to minimize average cost, with the 
previously hard (deterministic) constraints becoming chance constraints as explained next. 



Redefining the line flow constraints 

Power lines do not fail (i.e., trip) instantly when their flow limits are exceeded. A line carrying flow that exceeds 
the line's thermal limit will gradually heat up and possibly sag, increasing the probability of an arc (short circuit) 
or even a contact with neighboring lines, with ground, with vegetation or some other object. Each of these events 
will result in a trip. The precise process is extremely dificult to callibrate (this would require, among other factors, 
an accurate representation of wind strength and direction in the proximity of the line) [50| . Additionally, the rate at 
which a line overheats depends on its overload which may dynamically change (or even temporarily disappear) as 
flows adjust due to external factors; in our case fluctuations in renewable outputs. What can be stated with certainty 
is that the longer a line stays overheated, the higher the probability that it will trip - to put it differently, if a line 
remains overheated long enough, then, after a span possibly measured in minutes, it will trip. In summary, (thermal) 
tripping of a line is primarily governed by the historical pattern of the overloads experienced by the line, and thus it 
may be influenced by the status of other lines (implicitly, via changes in power flows); further, exogenous factors can 
augment the impact of overloads. 

Even though an exact representation of line tripping seems difficult, we can however state a practicable alternative. 
Ideally, we would use a constraint of the form "for each line, the fraction of the time that it exceeds its limit within a 
certain time window is small" . Direct implementation of this constraint would require resolving dynamics of the grid 
over the generator dispatch time window of interest. To avoid this complication, we propose instead the following 
static proxy of this ideal model, a chance constraint: we will require that the probability that a given line will exceed 
its limit is small. Let fij be the flow on line {i, j}, where the bold face indicates that it is a random quantity. Denote 
the "small probability" above by ey, and the flow limit on line {i, j} by f*™ ax ■ Then the chance constraint for each 
line, {i, j}, is: 

CCTL: V P(\f tJ \ >m ax ) <e ir (6) 
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Chance constraints [iff, [17|, [18| are but one possible methodology for handling uncertain data in optimization. 
Broadly speaking, this methodology fits within the general field of stochastic optimization. Constraint (j6)) can be 
viewed as a "value-at-risk" statement; the closely-related "conditional value at risk" concept provides a (convex) 
alternative, which roughly stated constrains the expected overload of a line to remain small, conditional on there 
being an overload (see [9| for definitions and details). 

Yet another alternative would be to rely on robust optimization. In this setting we would view the output of a 
renewable source j as (for example) an unknown quantity constrained to lie in an interval [lj,Uj], and formulate an 
optimization problem that requires (for example that the flow of each line stay within its limit no matter what the 
value of each renewable output j is, so long as it remains in its corresponding interval [lj,Uj]. Other variations are 
possible, but always with the same deterministic flavor. See e.g. This general approach is attractive because it 

makes few assumptions on the underlying uncertain process. However, we would argue that our chance constrained 
approach is reasonable (in fact: compelling) in view of the nature of the line tripping process we discussed above. 



Uncertain power sources 



We assume that there is a collection of wind sources (farms) , with one wind source located at each node in a given 
subset W of V. For each i £ W, the amount of power generated by source i is assumed to be of the form /x, + Wi, 
where /Ltj is constant, assumed known from the forecast, and Wi is a zero mean independent random variable with 
standard deviation oi. 

The physical assumptions behind this model of uncertainty are as follows. Independence of fluctuations at different 
sites is due to the fact that the wind farms are sufficiently far away from each other. For the typical OPF time span of 
15 min and typical wind speed of lOm/s, fluctuations of wind at the farms more than lOfcm apart are not correlated. 
We also rely on the assumption that transformations from wind to power at different wind farms is not correlated. 

To formulate and calibrate our models, we make simplifying assumptions that are approximately consistent with 
our general physics understanding of fluctuations in atmospheric turbulence; in particular we assume Gaussianity of 
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511 ] . We will also assume that only a standard weather forecast (coarse-grained on minutes and kilometers) is 



available, and no systematic spillage of wind in its transformation to power is applied 52]. 

We also have a strong additional -and purely computational- reason for the Gaussian assumption. As we will 
show below, under this assumption chance constraint ^ can be captured using a tractable deterministic optimization 
problem; our substitute for standard OPF (see below for details). As we will indicate below, in the case of non- 
Gaussian distributions our approach is easily modified so as to retain tractability, but at the cost of relying on a 
conservative approximation of the chance constraint. 

Other fitting distributions considered in the wind-modeling literature, e.g. Weibull distributions and logistic distri- 
butions 2^, 21 1, will be discussed later in the text as well [53j ] . In particular, we will demonstrate on out-of-sample 
tests that the computationally advantageous Gaussian modeling of uncertainty allows as well to model effects of other 
distributions. Our approach relies on a "robust" wind forecast; typically this would involve obtaining a reasonably 
accurate estimate on mean wind strength and a conservative variance estimation at each farm. This robust forecast 
will be used to compute our control as indicated in the next section. The overall edifice should be validated with 
actual data or out-of-sample simulations. (See Discussion and Methods for details.) 



Affine Control 



Since the power injections at each bus are fluctuating, we need a control to ensure that generation is equal to 
demand at all times within the time interval between two consecutive OPFs. We term the joint result of the primary 
control and secondary control the affine control. The term will intrinsically assume that all governors involved in 
the controls respond to fluctuations in the generalized load (actual demand which is assumed frozen minus stochastic 
wind resources) in a proportional way, however with possibly different proportionality coefficients: 

Vnodei e S : Pi = Pi - OH ^ u)j. cti > 0, (7) 

jew 

Here the quantities pi > and ai > are design variables satisfying (among other constraitns) YlieS a i = 1- Notice 
that we do not set any m to a standard (fixed) value, but instead leave the optimization to decide the optimal value. 
(In some cases it may even be advantageous to allow negative a, but we decided not to consider such a drastic change 
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of current policy in this study.) The generator output pi combines a fixed term pi and a term which varies with wind, 
—cti J2jew ■ Observe that J^i Pi = J2i Pi ~ Si w i s that is, the total power generated equals the average production 
of the generators minus any additional wind power above the average case. 

This affine control scheme creates the possibility of requiring a generator to produce power beyond its limits. With 
unbounded wind, this possibility is inevitable, though we can restrict it to occur with arbitrarily small probability, 
which we will do with additional chance constraints for all controllable generators, Vg G 9, 

CCgen: P(pf " < p g - a g £ Wj < p™*) > 1 - e g . (8) 

jew 



CC-OPF: Formal Expression 

We can now formally state the main optimization problem we introduce in this manuscript. CC-OPF is given by 
the following modification of the standard OPF problem ([I} : 

CC-OPF: minE w {c(jp)] s.t. (9) 

PFAV: Power Flow Eqs. (|2|) under average wind 
CCTL: Chance constrained line limits, Eq. © 
CCgen: Chance constrained generator bounds, Eq. © 

where E w [c(p, a)] is the expected cost of generation over the varying wind power w. If the cost function c is convex 
quadratic (standard practice), then so is this expectation. The expected cost objective and the chance constrained 
condition in CC-OPF both account for fluctuations in wind and also for the standard generation adjusting to these 
fluctuations via the aforementioned proportional control. 

p2j considers the standard OPF problem under stochastic demands, and describes a method that computes fixed 
generator output levels to be used throughout the period of interest, independent of demand levels. In order to handle 
variations in demand, instead relies on the concept of a slack bus. A slack bus is a fixed node that is assumed 
to compensate for all generation/demand mismatches - when demand exceeds generation the slack bus injects the 
shortfall, and when demand is smaller than generation the slack bus absorbs the generation excess. A vector of 
generations is acceptable if the probability that each system component operates within acceptable bounds is high 
- this is a chance constraint. To tackle this problem [22j proposes a simulation-based local optimization system 
consisting of an outer loop used to assess the validity of a control (and estimate its gradient) together with an inner 
loop that seeks to improve the control. Experiments are presented using a 5-bus and a 30-bus example. 

Chance constrained optimization has also been discussed recently in the Unit Commitment setting - discrete-time 
planning for operation of large generation units on the scale of hours-to- months accounting for the long-term wind- farm 
generation uncertainty [23l425l |. 

At first glance CC-OPF constitutes a significantly more difficult problem than the original OPF. Not only does this 
formulation contain additional (affine control) variables not present in the standard OPF, but even more significantly 
the chance constraints and the objective render it into a stochastic optimization problem, requiring the evaluation of 
expectation and the probability of overload over the exogenous statistics of the wind. A direct computational approach 
to solving problem ©, as in [22[, though universal (applicable to any type of exogenous distribution), would require a 
number of technical assumptions and elaborations to guarantee convergence and feasibility and would be prohibitively 
expensive, e.g. as discussed in (2^ |. Our technical approach is given next. 



CC-OPF: From Stochastic Formulation to Conic Programming 



Our methodology applies and develops general ideas of [9j to the power engineering setting of the generation re- 
dispatch under uncertainty. We show that under the assumptions of the basic power flow linearity, proportionality 
of the standard generation response to fluctuations, along with Gaussianity and statistical independence of the wind 
fluctuations at different sites, the CC-OPF © is reduced exactly to a deterministic optimization problem, see Eq. (flU)) 
of the Methods. Moreover this deterministic optimization over p and a is a convex optimization problem, more pre- 
cisely, a Second-Order Cone Program (SOCP) [IMIi^j allowing an efficient computational implementation discussed 
in detail in the "Cutting Plane" Subsection of the Appendix. 



8 



Let us emphasize that many of our assumptions leading to the computational efficient formulation are not restrictive 
and allow natural generalizations. In particular, using techniques from it is possible to relax the phenomenolog- 
ically reasonable but approximately validated assumption of wind source Gaussianity (validated according to actual 
measurements of wind, see pol |2~0] and references therein). For example, using only the mean and variance of output 
at each wind farm, one can use Chebyshev's inequality to obtain a similar though more conservative formulation. And 
following Q we can also obtain convex approximations to (J5|) which are tighter than Chebyshev's inequality, for a 
large number of empirical distributions discussed in the literature. In any case, we will perform (below) out-of-sample 
experiments involving our controls; first to investigate the effect of parameter estimation errors in the Gaussian case, 
and, second, to gauge the impact of non-Gaussian wind distributions. 

EXPERIMENTS/RESULTS 

Here we will describe qualitative aspects of our affine control on small systems; in particular we focus on the contrast 
between standard OPF and CC-OPF, on problematic features that can arise because of fluctuating wind sources and 
on out-of-sample testing of the CC-OPF solution, including the analysis of non-gaussian distributions. Some of our 
tests involve the BPA grid, which is large; in Supplementary Information we also present a second set of tests that 
address the scalability of our solution methodology to large cases. 

Failure of standard OPF 

Above (see cq. (fT])) we introduced the so-called standard OPF method for setting traditional generator output 
levels. When renewables are present, the natural extension of this approach would make use of some fixed estimate 
of each rencwablc's output (e.g., its mean output) and to handle fluctuations in renewable output through the same 
method used to handle changes in load: ramping output of traditional generators up or down in proportion to the net 
increase or decrease in renewable output. This feature could seamlessly be handled using today's control structure, 
with each generator's output adjusted at a fixed (preset) rate. For the sake of simplicity, in the experiments below 
we assume that all ramping rates are equal. 

Different assumptions on these fixed rates will likely produce different numerical results; however, this general 
approach entails an inherent weakness. The key point here is that mean generator output levels as well as in 
particular the ramping rates would be chosen without considering the stochastic nature of the renewable output 
levels. Our experiments are designed to highlight the limitations of this "risk-unaware" approach. In contrast, our 
CC-OPF produces control parameters (the p and the a) that are risk-aware and, implicitly, also topology-aware - in 
the sense of network proximity to wind farms. 

We first consider the IEEE 118-bus model with a quadratic cost function, and four sources of wind power added 
at arbitrary buses to meet 5% of demand in the case of average wind. The standard OPF solution is safely within 
the thermal capacity limits for all lines in the system. Then we account for fluctuations in wind assuming Gaussian 
and site-independent fluctuations with standard deviations set to 30% of the respective means. The results, which 
are shown in Fig. [21 illustrate that under standard OPF five lines (marked in red) frequently become overloaded, 
exceeding their limits 8% or more of the time. This situation translates into an unacceptably high risk of failure for 
any of the five red lines. This problem occurs for grids of all sizes; in Fig. [3] we show similar results on a 2746-bus 
Polish grid. In this case, after scaling up all loads by 10% to simulate a more highly stressed system, we added wind 
power to ten buses for a total of 2% penetration. The standard solution results in six lines exceeding their limits over 
45% of the time, and in one line over 10% of the time. 

Cost of reliability under high wind penetration 

The New York Times article "Wind Energy Bumps Into Power Grid's Limits" [28[ discusses how transmission line 
congestion has forced temporary shutdown of wind farms even during times of high wind. Our methodology suggests, 
as an alternative solution to curtailment of wind power, an appropriate reconfiguration of standard generators. If 
successful, this solution can use the available wind power without curtailment, and thus result in cheaper operating 
costs. 

As a (crude) proxy for curtailment, we perform the following experiment, which considers different levels of renewable 
penetration. Here, the mean power outputs of the wind sources are kept in a fixed proportion to one another and 
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FIG. 2: 118-bus case with four wind farms (green dots; brown are generators, black are loads). Shown is the standard OPF 
solution against the average wind case with penetration of 5%. Standard deviations of the wind are set to 30% of the respective 
average cases. Lines in red exceed their limit 8% or more of the time. 




FIG. 3: Failure of the standard OPF shown for a snapshot of the standard OPF solution on Polish grid from MATPOWER; 
full snapshot is shown in Section III of SI. Color coding and conditions of the experiment are equivalent to these of Fig. [2] 

proportionally scaled so as to vary total amount of penetration, and likewise with the standard deviations. First, 
we run our CC-OPF under a high penetration level. Second, we add a 10% buffer to the line limits and reduce 
wind penetration (i.e., curtail) so that under the standard OPF solution line overloads are reduced to an acceptable 
level. Assuming zero cost for wind power, the difference in cost for the high-penetration CC-OPF solution and the 
low-penetration standard solution are the savings produced by our model (generously, given the buffers). 

For the 39-bus case, our CC-OPF solution is feasible under 30% of wind penetration, but the standard solution 
has 5 lines with excessive overloads, even when solved with the 10% buffer. Reducing the penetration to 5% relieves 
the lines, but more than quadruples(l) the cost over the CC-OPF solution. See Figure |U Note that this approach 
does not only show the advantage of the CC-OPF over standard OPF but also provides a quantitative measure of the 
advantage, thus placing a well-defined price tag on reliability. 
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FIG. 4: 39-bus case under standard solution. Even with a 10% buffer on the line flow limits, five lines exceed their limit over 
5% of the time with 30% penetration (left). The penetration must be decreased to 5% before the lines are relieved, but at great 
cost (right). The CC-OPF model is feasible for 30% penetration at a cost of 264,000. The standard solution at 5% penetration 
costs 1,275,020 - almost 5 times as much. 



FIG. 5: 39-bus case. Red lines indicate high probability of flow exceeding the limit under the standard OPF solution. Generators 
are shades of blue, with darker shades indicating greater absolute difference between the chance-constrained solution and the 
standard solution. 



We have established that under fluctuating power generation, some lines may exceed their flow limits an unaccept- 
able fraction of the time. Is there a simple solution to this problem, for instance, by carefully adjusting (a posteriori 
of the standard OPF) the outputs of the generators near the violated lines? The answer is no. Power systems exhibit 
significant non-local behavior. Consider Fig. [5] In this example, the major differences in generator outputs between 
the standard OPF solution and our CC-OPF model's solution are not obviously associated with the different line vio- 
lations. In general, it seems that it would be difficult to by-pass CC-OPF and make small local adjustments to relieve 
the stressed lines. On the positive side, even though CC-OPF is not local and requires a centralized computation, it 
is also only slightly more difficult than the standard OPF in terms of implementation. 




Non-locality 
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FIG. 6: 39-bus case with four wind farms (green dots; brown are generators, black are loads). Lines in red are at the maximum 
of e = .02 chance of exceeding their limit. The three cases shown are left to right .1%, 8%, and 30% average wind penetration. 
With penetration beyond 30% the problem becomes infeasible. 




Increasing penetration 

Current planning for the power system in the United States calls for 30% of wind energy penetration by 2030 
Investments necessary to achieve this ambitious target may target both software (improving operations) and 
hardware (building new lines, sub-stations, etc), with the former obviously representing a much less expensive and 
thus economically attractive option. Our CC-OPF solution contributes toward this option. A natural question that 
arises concerns the maximum level of penetration one can safely achieve by upgrading from the standard OPF to our 
CC-OPF. 

To answer the question we consider the 39-bus New England system (from pfjl ]) case with four wind generators 
added, and line flow limits scaled by .7 to simulate a heavily loaded system. The quadratic cost terms are set to 
rand(0,l) + .5. We fix the four wind generator average outputs in a ratio of 5/6/7/8 and standard deviations at 30% 
of the mean. We first solve our model using e = .02 for each line and assuming zero wind power, and then increase 
total wind output until the optimization problem becomes infeasible. See Figure [6l This experiment illustrates that 
at least for the model considered, the 30% of wind penetration with rather strict probabilistic guarantees enforced by 
our CC-OPF may be feasible, but in fact lies rather close to the dangerous threshold. To push penetration beyond 
the threshold is impossible without upgrading lines and investing in other (not related to wind farms themselves) 
hardware. 



Changing locations for wind farms 

In this example we consider the effect of changing locations of wind farms. Wc take the MATPOWER 30-bus 
case with line capacities scaled by .75 and add three wind farms with average power output in a ratio of 2/3/4 and 
standard deviations at 30% of the average. Two choices of locations are shown in figure[71 The first remains feasible for 
penetration up to 10% while the second can withstand up to 55% penetration. This experiment shows that choosing 
location of the wind farms is critical for achieving the ambitious goal of high renewable penetration. 
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FIG. 8: 9-bus case, 25% average penetration from two wind sources. With shifting winds, the flow on the orange line changes 
direction with a large absolute difference. 

Reversal of line flows 

Here we consider the 9-bus case with two wind sources and 25% average penetration and standard deviations set 
to 30% of the average case and analyze the following two somewhat rare but still admissible wind configurations: (1) 
wind source (a) produces its average amount of power and source (b) three standard deviations below average; (2) 
the reverse of the case (1). This results in a substantial reversal of flow on a particular line shown in Figure [5J This 
example suggests that when evaluating and planning for grids with high-penetration of rencwables one needs to be 
aware of potentially fast and significant structural rearrangements of power flows. Flow reversals and other qualitative 
changes of power flows, which are extremely rare in the grid of today, will become significantly much more frequent 
(typical) in the grid of tomorrow. 

Out-of-sample tests 

We now study the performance of our method when there are errors in the underlying distribution of wind power. 
We consider two types of errors: (1) the true distribution is non-Gaussian but our Gaussian fit is "close" in an 
appropriate sense, and (2) the true distribution is Gaussian but with different mean or standard deviation. The 
experiments in this section use as data set the BPA grid, which as noted before has 2209 buses and 2866 lines, and 
collected wind data; altogether constituting a realistic test-case. 

We first consider the non-Gaussian case, using the following probability distributions, all with fatter tails than 
Gaussian: (1) Laplace, (2) logistic, (3) Wcibull (three different shapes), (4) t location-scale with 2.5 degrees of 
freedom, (5) Cauchy. For the Laplace and logistic distributions, we simply match the mean and standard deviation. 
For the Wcibull distribution, we consider shape parameters k = 1.2,2,4 and choose the scale parameter to match the 
standard deviations. We then translate to match means. For the t distribution, we fix 2.5 degrees of freedom and 
then choose the location and scale to match mean and standard deviation. For the Cauchy distribution, we set the 
location parameter to the mean and choose the scale parameter so as to match the 95th percentiles. 

We use our model and solve under the Gaussian assumption, seeking a solution which results in no line violations 
for cases within two standard deviations of the mean, i.e. a maximum of about 2.27% chance of exceeding the 
limit. We then perform Monte Carlo tests drawing from the above distributions to determine the actual rates of 
violation. See Figure [9] The worst-performer is the highly-asymmetric (and perhaps unreasonable) Wcibull with 
shape parameter 1.2, which approximately doubles the desired maximum chance of overload. Somewhat surprisingly, 
the fat-tailed logistic and Student's t distributions result in a maximum chance of overload significantly less than 
desired, suggesting that our model is too conservative in these cases. 

Next we consider the Gaussian case with errors. We solve with nominal values for the mean and standard deviation 
of wind power. We then consider the rate of violation after scaling all means and standard deviations (separately) . 
While the solution is sensitive to errors in the mean forecast, the sensitivity is well-behaved. With a desired safety 
level of e = 2.27% for each line, an error in the mean of 25% results in a maximum 15% chance of exceeding the limit. 
The solution is quite robust to errors in the standard deviation forecast, with a 25% error resulting in less than 6% 
chance of overload. See Figure [TUJ 
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Distribution Max. prob. violation 

Normal 0.0227 

Laplace 0.0297 

logistic 0.0132 

Weibull, k = 1.2 0.0457 

Weibull, k = 2 0.0355 

Weibull, k = 4 0.0216 

t location-scale, v = 2.5 0.0165 

Cauchy 0.0276 

FIG. 9: Maximum probability of overload for out-of-sample tests. These are a result of Monte Carlo testing with 10,000 
samples on the BPA case, solved under the Gaussian assumption and desired maximum chance of overload at 2.27%. 
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Forecast error 



FIG. 10: BPA case solved with average penetration at 8% and standard deviations set to 30% of mean. The maximum 
probability of line overload desired is 2.27%, which is achieved with forecast error on the graph. Actual wind power means 
are then scaled according to the x-axis and maximum probability of line overload is recalculated (blue). The same is then done 
for standard deviations (green). 



DISCUSSIONS 



This manuscript suggests a new approach to incorporating uncertainty in the standard OPF setting routinely used 
in the power industry to set generation during a time window, or period (typically 15 min to one hour duration). 
When uncertainty associated with renewable generation is quantified in terms of the probability distribution of output 
during the next period, we incorporate it through chance constraints - probabilistic conditions which require that any 
line of the system will not be overloaded for all but a small fraction of time (at most one minute per hour, for 
example). Additionally, the modeling accounts for local frequency response of controllable generators to renewable 
changes. The key technical result of this manuscript is that the resulting optimization problem, CC-OPF, can be 
stated as a convex, deterministic optimization problem. This result also relies on plausible assumptions regarding the 
exogenous uncertainty and linearity of the underlying power flow approximations/equations. In fact, our CC-OPF is 
a convex (conic) optimization problem, which we solve very efficiently, even on realistic large-scale instances, using a 
sequential linear cutting plane algorithm. 

This efficient CC-OPF algorithm becomes an instrument of our (numerical) experiments which were performed on 
a number of standard (and nonstandard) power grid data sets. Our experimental results are summarized as follows: 

• We observe that CC-OPF delivers feasible results where standard OPF, run for the average forecast, would fail 
in the sense that many lines would be overloaded an unacceptably large portion of time. 

• Not only is CC-OPF safer than standard OPF, but it also results in cheaper operation. This is demonstrated by 
considering the optimal cost of CC-OPF under sufficiently high wind penetration solution (where standard OPF 
would fail) and the low penetration solution (corresponding to the highest possible penetration where standard 
OPF would not fail). 

• We discover that solutions produced by CC-OPF deviate significantly from what amounts to a the naive adjust- 
ment of the standard OPF obtained by correcting dispatch just at those generators which are close to overloaded 
lines. 

• We test the level of wind penetration which can be tolerated without upgrading lines. This experiment illustrates 
that, at least for the model tested, the 30% of wind penetration with rather strict probabilistic guarantees 
enforced by our CC-OPF may be feasible; but much lower wind penetration remains feasible under the standard 
approach. 
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• We experiment with location of wind farms and discover strong sensitivity of the maximum level of penetration 
on the choice of location - optimal choice of wind farm location is critical for achieving the ambitious goal of 
high renewable penetration. 

• Analyzing fluctuations of line flows within CC-OPF solution admissible under high wind penetration, we discover 
that these fluctuations may be significant, in particular resulting in reversal of power flows over some of the lines. 
This observation suggests that flow reversals and other qualitative changes of power flows, which are extremely 
rare in the grid of today, will become significantly much more frequent (typical) in the grid of tomorrow. 

• We also studied an out-of-sample test consisted in applying CC-OPF (modeling exogenous fluctuations as 
Gaussian) to other distributions. Overall these tests suggest that with a proper calibration of the effective 
Gaussian distribution our CC-OPF delivers a rather good performance. One finds that the worst CC-OPF 
performance is observed for the most asymmetric distributions. 

The nature of the problem discussed in the manuscript - principal design of a new paradigm for computationally 
efficient generation re-dispatch that accounts for wind fluctuations - inevitably required incorporation of a number of 
assumptions and approximations. In particular, we made simplifying assumptions about static forecasts and general 
validity of power flow linearization. We have also focused solely on failures associated with line congestion ignoring 
other possible difficulties, for example those associated with loss of synchronicity and voltage variations. However, 
we would like to emphasize that all of these assumptions made (admittedly natural for a first attack on the problem) 
also allow generalization within the approach just sketched: 

• Accounting for time evolving forecast /loads/etc. Wind forecast, expressed in terms of the mean and standard 
deviation at the wind farm sites, changes on the scales comparable to duration of the generation re-dispatch 
interval. Loads may also change at these time scales. When the slowly evolving, but still not constant, wind 
and load forecasts are available we may keep the quasi-static power flow description and incorporate this slow 
evolution in time into the chance constrained scheme. These changes will simply result in generalizing the 
conic formulation of Eq. [10] by splitting what used to be a single time interval into sub-intervals and allowing 
the regular generation to be re-dispatched and parallel coefficients to be adjusted more often. Ramping rate 
constraints on the controllable generation may naturally be accounted in the temporal optimization scheme as 
well. 

• Accounting for nonlinearity in power flows. Evolution of the base case invalidates the linearization (DC-style) 
hypothesis. However, if variations around one base case becomes significant one may simply adjust the lin- 
earization procedure doing it not once (as in the case considered in the manuscript) but as often as needed. 
Slow adjustment of the base case may also be included into the dynamical procedure mentioned one item above. 
Additionally, some interesting new methodologies for handling nonlinearities have recently emerged, see (30j . 

• Accounting for synchronization bounds. Loss of synchronicity and resulting disintegration of the grid is probably 
the most acute contingency which can possibly take place in a power system. The prediction of those conditions 
under which the pow er g rid will lose synchronicity is a difficult nonlinear and dynamic problem. However, 
as shown recently in [3lJ, one can formulate an accurate linear and static necessary condition for the loss 
of synchronicity. A chance-constrained version of the linear synchronization conditions can be incorporated 
seamlessly in our CC-OPF framework. 

Finally, we see many opportunities in utilizing the CC-OPF (possibly modified) as an elementary unit or an integral 
part of even more complex problems, such as combined unit commitment (scheduling large power plants normally 
days, weeks or even months ahead) [HI with CC-OPF, planning grid expansion [32| while accounting for cost operation 



under uncertainty, or incorporating CC-OPF in mitigating emergency of possible cascades of outages [33H40J. In this 



context, it would be advantageous to speed up our already very efficient CC-OPF even further. See, for example, [41 



42] . A different methodology, relying on distributed algorithms, can be found in |4S 



METHODS 

CC-OPF as Deterministic (Conic) Programming 



Eq. (j9]) states the problem of generation dispatch under uncertainty due to wind as a stochastic optimization problem 
with chance (probabilistic) constraints. The critical part of our approach consists in providing a convex expressions for 
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the expectation of the objective and for the probabilities in Eq. ©, as a function of generation dispatch optimization 
variables. This gives rise to the following deterministic optimization problem (see Section I of SI for derivations and 
proofs) 



mm 



S.t. 



ca Pi + a? °i + Cl2p * + Ca ( 10 ) 
ies \ \ jew J J 

^ / a i = l, a>0, p>0 (11) 

B9 = p + [i-d, Y,(Pi + W - *) = (12) 

iev 

ri-1 

for 1 < i < n - 1 : ^ BySj = a,, S n = 0, (13) 
j'=i 

V{i,i}G£: 

4 > 4 E - B tk S t + 5rf, (14) 

\Pii(fii M < f™* - <j>~\l - tij/2) Sij (15) 

The objective here is simply K w [c(p, a)] written explicitly; the vectors p and a model our control methodology as 
described above; 9 is the average phase angle vector; the first equation in (|12|) amounts to the standard DC model 
equation relating phase angle to power injections. The second constraint in [12] is a flow balance statement: the sum 
of all power injections is zero. The nature of our control will guarantee that power is balanced for any configuration 
of wind power. 

Constraints (jT4j) and (|15|) constitute a (deterministic) representation of the chance constraint (|6]). Here, /% is 
the susceptance on line <Jh is the standard deviation of wind source k; </> is the standard normal cumulative 

distribution function; B + is an appropriate, sparse generalized inverse of the bus susceptance matrix; the 5 are 
auxiliary variables. It can be shown (see SI) that, under the indepencence assumption for the wind fluctuations Wi, 
but without requiring Gaussianity, the right-hand side of constraint (1141) is the variance of the flow on line (i, j). Thus, 
at optimality (p~4|) and (fTS"]) will amount to 

< fT x - ^C 1 - W 2 K' ( 16 ) 

where <7y- is the standard deviation of the flow fij on line («,j). 

Further, it can be shown that the bus angles Qi and line flows fij both are affine combinations of the Wi. Thus, 
assuming that the Wi arc Gaussian, then so will be the 9i and fij. Since the left-hand side of equation (fT5|) is the 
absolute value of the expected flow on line it follows under the Gaussianity assumption that (|16[) is as claimed 

a valid representation of the chance constraint ([6]): it states that the expectation of flow on line is the right 

multiple of a standard deviation away from the maximum f^ ax , as per the risk tolerance . 

In deriving (|15j) we are thus explictly making use of the Gaussianity assumption. However, using the results in Q one 
can show that in the case of arbitrary distributions of the Wi (but assuming finite variances) one can obtain a convex 
(conic) conservative approximation to the chance constraint by simply replacing the quantity _1 (1 — £ij/2) in (|15[) 
with an appropriate coefficient f2 (effectively, this approach relies on estimates from the theory of large deviations). 
Thus, even in the non-Gaussian case, our general approach remains essentially the same. 



Cutting-Plane Algorithm 

The number of conic constraints (|14[) is equal to the number of lines, and in the case of a large grid this can 



prove challenging. For example, in the Polish 2003-2004 winter peak case, [44[ reports over 6000 conic constraints 
after pre-solving. All of the commercial solvers (4414461 we experimented with reported numerical difficulties and were 
unable to solve a Second Order Conic Programming (SOCP) of this size. 

In order to solve these large cases, we employed a cutting-plane algorithm. Refer to the Section 2 of SI for a detailed 
description. The algorithm repeatedly solves linear programs that include all the linear constraints in our formulation 
together with a finite set of first-order approximations to the conic constraints. Having solved one such problem the 
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algorithm checks for a violated conic constraint and if found, it approximates that conic constraint with its first-order 
approximation at the point of violation. We repeat this process until the largest violation of a conic constraint is 
sufficiently small, in our case less than 10 -6 . This method is able to solve cases with thousands of buses usually in 
just CPU seconds on a laptop computer. The following table shows typical behavior: 



Iteration Max rel. error 


Objective 


1 


1.2c-l 


7.0933e6 


4 


1.3e-3 


7.0934e6 


7 


1.9e-3 


7.0934e6 


10 


1.0e-4 


7.0964e6 


12 


8.9e-7 


7.0965e6 



Each row of this table shows that maximum relative error and objective value at the end of several iterations. The 
total runtime was 25 seconds. Note the "flatness" of the objective. This makes the problem nontrivial - the challenge 
is to find a feasible solution (with respect to the chance constraints); at the onset of the algorithm the computed 
solution is quite infcasiblc and it is this attribute that is quickly improved by the cutting-plane algorithm. 



Power Grid Models 

The BPA grid data was extracted from [13] (available to the authors) . Corresponding wind data was extracted from 
[l|. All other power grid cases, in particular Polish Grid, 118 bus, 39 bus and 30 bus systems are publicly available 
from [29I ]. The wind-related modifications for each example are explained in the text above. 

We are thankful to S. Backhaus and R. Bent for their comments, and to K. Dvijotham for help with BPA data. 
The work at LANL was carried out under the auspices of the National Nuclear Security Administration of the U.S. 
Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396. MC and SH 
also acknowledge partial support of the DTRA Basic Research grant BRCALL08-Pcr3-D-2-0022. DB was partially 
supported by DOE award Ide : sc0002676l 



Supplementary Information: CC-OPF as Conic Programming 

Notation 

The following notation and conventions will be used throughout, (i) n = number of buses, m = number of lines, 
(ii) Let B be the (n — 1) X (n — l)-matrix obtained by removing from B row and column n. (Assuming the grid is 
connected, B is invertible.) (iii) In what follows we use bold font to indicate uncertain quantities, (iv) We will write 
S = J2iew w i> a2 = Sigw °f ■ ( v ) F° r convenience of notation, here we assume that bus n does not hold a generator 
and does not hold a wind farm, in other words, n ^ 9 U W. This assumption is easily attained by adding a "dummy" 
bus (with zero generation, demand and wind output) and attaching it to the grid with a dummy line, (vi) For a bus 
i, let di be its demand. When i ^ D we write di = 0. Write pi = 014 = for each bus i that is not under the affinc 
control discipline. For each i with i £ W let ^ = and let Wi be the random variable with mean and variance equal 
to 0. Then, bi = bi + pi — aiS + Wi, where bi = fii — di, is the net power injected into the network at bus i. Note also 
that we must always have 

= ^ bi = y^(gj - QjjS) + ^2((M +Wi) - J^dj 



Analysis 

Let b = (bi, . . . ,&„_i) T , a = (an,. . .,a„_i) T , p = (pi,. . . ,p„_i) T , and w = (/ii + w 1} . ..fi n -i + u>n-i) T - Thus, 
given a choice of the control variables pi and on, we have that the (random) phase angle vector 9 = (0±, . . . , 9 n —i) T 
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satisfies: 

B6 = b + p - Sa + w. (18) 
The vector 6 has one degree of freedom, and so we set 9 n = and thus E(9 n ) = 0. So for all 1 < i < n, 

n-i {[B- l (Jb + p)]i, i < n, 

6 i = a i -5 l S+J2^w j , 6i=l (19) 
i =1 |0, otherwise. 

^ ^H-ET 1 ]^, i < n, s ^[[B- l a]i,i < n, ^ 
0, otherwise. * 1 0, otherwise. 



Consider a line (i, j). Since the flow fij on line equals Pij(0i — Oj) we have that 

n-l 

fij = Oi - 9j - (Si - Sj)S + ^(7r lfe - n jk )w k , (21) 



fe=i 



and therefore, E(fij) = (3ij(9i — 9j) and, since the tOi are pairwise independent, 



vartfa ) = 4 J2 "i(*ik - Kjk - Si + S,) 2 . (22) 
fcew 

Likewise, denoting by P g the power produced by generator g, we have E(P g ) = p g , var(P g ) = OigYljew a j ■ As a 
result of the above we have: 

Lemma .1 For a given choice of vectors p and a, each quantity 6i or fij is an affine function of the random 
variables Wi . 

Proof. Follows from eqs. (JTOJl and ([2"0J). ■ 

Formulation. Let the vectors p and a (both in R" _1 ) be fixed, and consider the following system of inequalities, on 
variables Si,9i (all for 1 < i < n — 1) and /i,j,Sj,j (for each line (i,j)) [we also have the additional quantities 5 n and 
9 n as variables fixed at zero, for convenience of notation]: 

n-i 

forl<t<n-l, ^ Bij Sj = an; S n = 0, (23) 
i=i 

V (i, j) : 4- ^ ^ (7r ife - 7T ifc - + <^) 2 < 4 , (24) 

for 1 < i < n - 1, ^3 ~Pi = h h = 0, (25) 

j'=i 

V(i,j): - - 0,) = 0. (26) 



Theorem .2 Consider the affine control given by a pair of vectors p and a satisfying Ji7| ). Then (6, 9, f, s) is feasible 
for ^3\) - H26\) if and only if: (a) for each bus i, 9i = E(9i); (b) for each line (i,j), fi.j = E(fi^); (c) for each line 
(i,j), sfj > var(f itj ); (d) for each generator g, E(P g ) = p g and var(P g ) = a 2 g J2jew °j • 

Proof. Suppose (5,6, f,s) is feasible for Eqs. dHSJ) - (EHJ) - By Eq. (g5]) we have that at 1 < i < n-l, 6~i = [B- 1 (b+p)] l = 
E(8i), and by Eqs. ©, (25), f itj = E(f i , j ), for each line (i,j). Similarly, Eq. (20]) , (H , ([22]) imply sfj > var(f i>:j ) 
as desired, (d) Holds by construction. The converse is similar. ■ 

Note: it is easily seen that under the conventions for bus n, constraints (11-14) of the main text are equivalent to 
Eqs. d23j)-((26l),(ITZ|)- 

We let < Eij < 1 (and, < e g < 1 ) denote the probabilistic tolerance for line (i, j) (resp., for generator g), and 
denote by e the (m + |S|)-vector of tolerances. 
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Definition .3 Given a vector e of tolerances, a control (p,a) is (1 — e)- strong, if Prob(\fij\ > f™ ax ) < tij, for each 
line (i,j)- and Prob(p™ ln < P g < p™ ax ^ > 1 — e g , for each generator g. 

For the next result, suppose we have a fixed control (p, a), and consider Eq. (fT9|) . By expanding S as J2k Wk > we 
see that for any line (i, j), fij equals a constant plus a linear combination of the Wk- It follows that if the Wk are 
normally distributed, then so is each flow value fij. 
Additional notation: For real < r < 1 we write g(r) = (/> _1 (1 — r). 

Lemma .4 Let (p,a) satisfy Eq. |-?7| ). Suppose (p,a) is (1 — e)-strong. Then there exists a vector (5,9, f,s) feasible 
for Eqs. 12:^ 1 - W(Ki such that, in addition, for all lines (i,j), 

i/«i</r x - 7 ?( e «) s «. ( 27 ) 

and for all generators g, 

pf n +v(e g )a g *]<p g <p™ aX -V(e g )a g ^. (28) 
Y jew Y jew 

Conversely, if there exists a vector (5,9, f,s) satisfying Eqs. 123\) -K2b)) and Eqs. |^7| j-fM|) then (p. a) is (1 — 2e)- strong. 

Proof. If Prob(\fij \ > f^ ax ) < e then using Theorem [2] we have that Eqs. (|27[) hold. For the converse, we have that 
Eq. (27]) implies Prob(fi,j > f% ax ) < e and Prob(f i:j < -f% ax ) < e. Likewise with Eqs. (28]). ■ 
Comment: Eqs. (|2"T]) . together with Eq. (2"o]) . are equivalent to the constraint Eq. (15) of the main text. 



Supplementary Information: Cutting-Plane Algorithm 

In the case of a grid with thousands of buses and lines, the optimization problem given by Eqs. ([IT]) , (j2lI]) - (f2l)]) . (j2"T]) - 
(|28[) (or, equivalently, Eqs. (10-15) of the main text), amounts to a large-scale convex conic programming problem. 
Experience with realistic examples with thousands of lines shows that commercial optimization packages are unable 
to solve the resulting problems "out of the box." Here we outline a simple algorithm that proves effective and fast, a 
so-called "cutting-plane" algorithm [26| . (4^ . From a theoretical standpoint the algorithm is motivated and justified 
by the "efficient separation is equivalent to efficient optimization" paradigm that underlies the ellipsoid method [48| . 

Without constraints (|24[) . (|28p . the optimization problem is a linearly constrained convex quadratic programming 
problem, still somewhat large (in the case of large grids) but within the reach of commercial solvers. Our algorithm 
iteratively replaces these constraints with a number of linear approximations which are algorithmically discovered. 

For a given line (i,j) write Cij(5) = Pij ^J2kew <T fc( 7r ifc — n jk ~ + ^j) 2 ! then constraint (24]) can be written as 
Cij (5) < Sij . Given a vector 5* , this constraint can be relaxed by the (outer) approximation CV,- (5* ) + dC Q$ S (5i — 
$i) + 9C dS S " — s ij> which is valid for all 5. A similar linearization applies to Eqs. ([2"5]) . This technique gives 

rise to the following iterative algorithm. The algorithm maintains a linear system of inequalities A(p, a, 5, 9, /, s) T > b, 
henceforth referred to as the master system, which is initialized to the set of all our constraints but with the conies 
removed; i.e. Eqs. (T7]),(I23]),(lS])-(in]),(E3- Denoting by F(p,a) the objective function in Eq. (10) of the main, the 
algorithm iterates through the following steps: 

1. Solve min{F(p, a) : A(p, a, 5, 9, f, s) T > b}. Let (p* ,a* ,5* ,9* , f* , s*) be an optimal solution. 

2. If all conic constraints are satisfied up to numerical tolerance by (p* , a* ,5* ,9* , /*, s*), Stop. 

3. Otherwise, add the outer inequality arising from that constraint (|24[) which is most violated to the master 
system. 

As the algorithm iterates the master system represents a valid relaxation of the constraints Eqs. (10-15) of the main 
text; thus the objective value of the solution computed in Step 1 is a valid lower bound on the value of problem. Each 
problem solved in Step 1 is a linearly constrained, convex quadratic program. Computational experiments involving 
large-scale realistic cases show that the algorithm is robust and rapidly converges to an optimum. 

Table |T] displays typical performance of the cutting-plane algorithm on (comparatively more difficult) large problem 
instances. In the Table, 'Polish2' is the case described in the main text (loads increased 10% as in the text). Polishl 
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and Polish3 are two others included in MATPOWER. All Polish cases have uniform random costs on [0.5, 2.0] for 
each generator and ten arbitrarily chosen wind sources. The average wind power penetration for the four cases is 
8.8%, 3.0%, 1.9%, and 1.5%. 'Iteratons' is the number of linearly-constrained subproblems solved before the algorithm 
converges. 'Barrier iterations' is the total number of iterations of the barrier algorithm in CPLEX over all subproblems, 
and 'Time' is the total (wallclock) time required by the algorithm. These instances all prove unsolvable if directly 
tackled by CPLEX or Gurobi. 

TABLE I: Performance of CC-OPF method on a number of large cases. 



Case 


Buses 


Lines 


Time (s) Iterations 


Barrier iterations 


BPA 


2209 


2866 


5.51 


2 


256 


Polishl 


2383 


2896 


13.64 


13 


535 


Polish2 


2746 


3514 


30.16 


25 


1431 


Polish3 


3120 


3693 


25.45 


23 


508 
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